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Abstract. In the past few y(^ars (:onsi<I(n'able progress has been made in Monte Carlo 
simulations of first-order phase t.raiisil.ioiis and in the analysis of the resulting finite-size 
data. In this pap(T sp(Mial (nnphasis will pla(:(Ml on multicauouical simiilal.ious using 
multigrid update techniques, on numerical estimates of interface tensions, and on accurate 
methods for determining the transition point and latent heat. 

1 Introduction 

First-order phase transitions play an important role in many fields of physics [1-3]. Well- 
known examples are field-driven transitions in magnets, crystal melting, the nematic- 

isotropic l.rausil.iou in lifiuid crysl.als or. at. much liiglior energy scales, tlio (locoiifiiiiiig 
transition in hot qiiark-gliioii matt(T and tli(^ various transitions in tli(^ (woliition of the 
early universe [.3]. 

An important property of first-order phase transitions is phase coexistence. For field- 
driven transitions as, e.g., in the Ising model at low temperatures or the (p''' theory dis- 
cussed below, this is refiected in the canonical ensemble by a highly double-peaked prob- 
ability distribution Pcan(™) of the magnetization m. To sample in Monte Carlo (MC) 
simuladous tlio two peaks willi (lie correct rolallvo weight llio system has to pass many 
tiuK^s through mixed phas(^ configurations. For finit(^ p(Tiodic systems of size L**, such con- 
figurations contain (at least) two interfaces whicli contribute an additional term 2aL'''^^ 
to the free energy, where a is the interface t(nisioii. Compared to the peak maxima they 
are hence suppressed by an additional Boltzmauu factor c< exp{— 2aL''~^} which implies 
an exponential divergence of the autocorrelation time with system size, r oc exp{2(TL''^^}. 
This property is sometimes called supercritical slowing down. The same arguments apply 
to temporatur("-(lriv(^i first-order phase transitions as, e.g. in Potts models, where m has 
to be repla(^ed by the energy c and Pcan(™) by Pa^nif)- 

Much effort has been spent in the past few years to develoj) officiout inothods for 
numerical studies of this important class of phase transitions. Both, improved update 
schemes for data generation and refined techniques for data analysis have been studied. 
To overcome the slowing down problem, a so-called multicanonical method has been 
proposed which is basically a reweighting approach and can, in principle, be combined 
with any legitimate update algorithm. To date most applications employed the local 
heat-bath or Metropolis algorithms. These studies showed that (at least for all practical 

*To appear in Computer Simulations in Condensed Matter Physics VII , eds. D.P. Landau, K.K. Mon, 
and H.B. Schiittler (Springer Verlag, Heidelberg, Berlin, 1994). 
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purposes) the exponential slowing down is indecMl nMluced to a mudi weaker power-law 
divergence with increasing system size. Since tlic remaining slowing down problem is, 
however, still sever(\ S(^v(n'al stii<li(^s hav(^ tried to combine the multicanonical approach 
with other update algorithms which are known to be much more efficient in the case of 
continuous phase transitions. In Sec. 2, after a brief review of the multicauouical molliod, 
special emphasis will be laid on a recently proposed combination with multigrid update 
algorithms which was shown to give a further real-time improvement of about one order 
of magnitude for a two-dimensional <f>* lattice model. 

The uniform accuracy of the probability distribution in multicanonical simulations led 
to r(iiiv(^stigatioiis of a long- known histogram technique to determin(^ th(^ iiit(^rfac(^ t(^ision 
between the coexisting phases at llio transition point. The 2D g-state Potts model, the 
2D and 3D Ising model and th(^ 2D (i* lattice model as well as the (3-|-l)D SU(3) lattice 
gauge theory have been studied so far. In Sec. 3.1 a summary is given of numerical results 
for the interface tension of 2D Potts models which can be compared with a recently 
derived analytical formula. This formula relies on an exact expression for the correlation 
length in the disordered phase at the transition point. Some direct numerical tests of this 
expression ai'(^ pr(V(^it(Ml in S(x:. 3.2. 

Parallel to the algorithmic developments many exact results coucoruiug the finite-size 
scaling b(^havior at (strong) first-order phase transitions could be prow^i. This formulation 
suggested refined methods to estimate the transition point and latent heat from finite- 
size data, which are discussed in Sec. 3.3. The important point is that these estimates 
are exponentially close to the infinite-volume limit, i.e., they do not show the typical 
power-law corrections oc 1/L'' of the traditional observables and are hence of improved 
accuracy. 

Finally, Sec. 4 contains a few concluding remarks. 

2 Improved Generation of Monte Carlo Data 
2.1 Multicanonical reweighting 

The slowing down problem at tirst-ord(T phas(5 transitions is din^'tly r(4at(Ml to the double- 
peak shape of the canonical probability distribution Pcan("j) or Pcan(e)- In multicanoni- 
cal simulations [4, 5] of field-driven transitions this problem is avoided by simulating an 
auxiliary distribution Pn„,cJm) = Pcanini) x exp(— /(?«,)), where the reweighting factor 
exp(— /(ni)) = w{rn)^^ is adjusted iteratively until Piii„<a("i) = const, between the two 
peaks [6-8]. This gives the mix(Hl phas(^ configurations tli(^ sam(^ statistical w(iglit as tli(^ 
pure phases. Precisely the same idea applies to temperature-driven transitions with in re- 
placed by ( . C'a.nonical (^xp(M'lalioii values (C')caii of any observable O can be recovered by 
the reweighting formula (C')can = ('i'C)muca/('w)niucai where (. . .)niuca denote expectation 
values with respect to the multicanonical distribution Pmuca- 

Using local algorithms to update the multicanonical distribution it was demonstrated 
for various models [9-16] that the exponential slowing down is indeed reduced to a power- 
law behavior r oc X''" = L''" with <i- ~ 1 — 1.5. as one would expect froiu a simple random 
walk argument. Wliih^ this is clearly an important st(^p forward. tli(^ r(unaiuing slowing 
down problem is still ficwrc. In fact, it is even wors(^ than in simulations of critical 
phenomena where standard local MC algorithms yield r oc with z ^ 2 [17, 18]. Here, 
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however, several improv(xl (mostly non-local) update aln'orilLnis (ov(^rr(^laxal.i()u. cluster, 
multigrid. . . . ) are available wliicli greatly reduce or even completely eliminate tlie slowing 
down pr()hl(nn [19-21]. 

In order to fiirtlier improve tlie performance of Monte Carlo sinnilations of first-order 
phase transitions it is therefore quite natural to study combinations of multicanonical 
reweighting with these improved update algorithms. Overrelaxation techniques have been 
used in the context of SU(3) lattice gauge theory [16], but due to the complexity of this 
system no systematic investigation of the performance has been reported. Cluster updates 
are difficult to implement in a straightforward way since, due to the reweighting factor 
f(m) or /(e), the multicanonical Hamiltonian is implicitly non-local. TIk^ more involved 
solution [11] proposed for Potts models yields a re<hice<l exjjonent a and should thus be 
favorable for large systems. For moderate lattic(^ siz(V. however, it is again not clear how 
much is gained in real computer time. Multigrid update techniques finally can be quite 
efficient for problems with continuous fields as will be discussed in the next subsections. 

2.2 Multigrid Monte Carlo Update Schemes 

The general strategy of multigrid MC update schemes [22-24] is to perform collective 
moves on successively long(n' l(nigtli scal(es in a systematic order as extensively discussed 
in the context of partial differential equations [25]. Both the type of collective moves 
and the sequence of length scales are parameters of multigrid schemes which can be 
defined in two equivalent ways: in a unigrid formulation where the update scheme is 
described in terms of the field variables on the original (fine) lattice, or in a recursive 
m.ultigrid formulation where successively coarser lattic(es and associated Hamiltonians are 
introduced [19]. While the unigrid formulation is very transparent and easy to program, 
the mor(e involwxl r(x:ursiv(e multigrid formulation is numerically much imm^ efhcient. 
In a unigrid formulation the collective update proposals are usually taken as square 

excitations, that is all fields in successively larger blocks of size 1, 2'', 4'' V = 2'"' are 

moved by the same amount <j>),. Another possibility are excitations of the form of pyramids 
where is the amplitude of the central displacement. We then accept or reject the 
proposed collective mov(> in accordance^ with the usual M(^.ropolis formula. In a recursive 
multigrid fornnilation the block amplitudes <j)), are identified with field variables ^f ' on 
coarse lattic(es S'*' (with 2*'' sit(es), and the shape of tlK^ (Excitations is controlled by an 
operator V which interpolates the coarse grid fields back to the next finer grid H'*"""^'. The 
squar(E (Excitations correspond to a piecewise constant interpolation ('P(0'*', 02*'i • ■ •) ~ 
((pf\(pf \ <j)2\<t>2 \ ■ ■ ■))> ™d the pyramids correspond to a piecewise linear interpolation. 
To update the fields (/i'*' one defines a coarse grid Hamiltonian recursively as i?'*' ((/)•*') — 
ij(t+i)(0(*+i) +-P(/,<''). iJl") = H. and initializ(Es th(! 4>f^ to zero. An update sweep 
over the coarse grid then corresponds to updating the amplitudes (f>h of blocks of 
size 2'"^'''' in the unigrid formulation. The multigrid formulation goi'ti actually one step 
further and also defines the sequence of coarse grids or length scales in a recursive way. 
More precisely the update on level k thus consist of a) mi presweeps using any valid 
update algorithm with Hamiltonian i?'*', b) calculating the Hamiltonian for the next 
coarser grid and initializing the field variables (jif^^^ to zero. One then c) updates the 
field variables (j)f^^'' by applying the multigritl update 7i_i times, d) interpolates the 
variables of level k — 1 back to the finer grid, and e) performs another m2 postsweeps 
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W cycle - a recursive "black box" 



W cycle (L=32) 




Figure 1: The recursive definition of the W-cycle and the expanded W-cycle for L = 32. 



on level k. On the coarsest grid S*"' only steps a) and e) can be performed and the 
recursion stops. At step c) the sequence of coarse grids is defined by the cycle parameters 
74;. For 74 = 7 = 1 (V-cycle) every grid is given the same weight, while for 74. = 7 = 2 
(W-cycle) the coarser grids are updated more frequently. The cycle names derive from 
their graphical representations which very much resemble the letters V and W; see Fig. 1. 

For the Ganssian mod(4 it can Ix^ prowni [19] that critical slowing down is c(mipl(Et(4y 
eliminated a) for piecewise linear interpolation (or higher) and any cycle or b) for piecewise 
constant interpolation and the W-cycle (or higher). For a V-cycle with piecewise constant 
interpolation this has not yet been proven and, in fact, numerical evidence indicates that 
in this case z = 1 [26, 27]. 

An important feature of the recursive multigrid formulation is the fact that the com- 
putational labor W per cycle is proportional to the work Wg per sweep of the original 
lattice. W !'l'o(l — 7/2'')^^ (7 < 2''). ind(5pendent of the lattice size. In the unigrid 
formulation, on the other hand, the estimates are W Ri Wg logj L for the V-cycle and 
W WoL'°i^^'' for 7 > 2, i.e, W WqL for a W-cycle. 

2.3 Multicanonical Multigrid Monte Carlo Method 

When adapting multigrid schemes to multicanonical simulations one only has to make 
sure that also the reweighting factor f(m) or /(e) is treated properly [27-29]. Let us 
first consider the case of a field-driven transition in the unigrid viewpoint. Using piece- 
wisc constant interpolation, a proposed move <f>i, for a block of size 2*"' would change the 
magn(Etization by Am — 2*'' (/>;,/!/. The decision of acceptance is thus simply to be based 
on the value of AE^„^^ = AiJ^an + /("'■ + '^'''''tpiV) — J(m). In the recursive nmltigrid 
implementation the multicanonical modification is very similar since the coarse grid mag- 
netization can be computed in the standard way. For a temperature-driven transition 
Am would have to be replaced by AE^^^JV . which is also straightforward to imph^iKEnt 
in both the unigrid and multigrid formulation. For a discussion of more general situations, 
see Ref.[29]. 
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2.4 Application to the 2D (f)* Model 



As an application of the multicanonical multigrid algorithm we considered in Ref. [29] the 
scalar 4>* lattice field theory defined by the partition function 

Z-U [/'^ exp (- E (i(V0.)' - + #^)) , (1) 

with g = 0.25 and /i^ > /(^((?) > 0. Here is a line of critical points S(^paratinn' tli(^ 

disordered {ft^ < fJ-Kg)) and ordered [fi^ > nKg)) phase. For d > 2 and fi^ > lil(g), the 
refiection symmetry is spontaneously broken in the infinite volume limit L — > oo, and the 
average magnetization (m) — ((1/y) J2i=i 4'i) acquires a nonvanishing expectation value. 
Consequently, if a term h (pi is added to the energy, the system exhibits first-order 
phase transitions driven by the field h. 

We have studiod the field-driven first-order phase trausitaou between the two ordered 
phases in d = 2 dinK^isions at // = 0.25 for fi^ = 1.. 30, 1.35, and 1.40. These points 
are sufficiently far away from ilic critical point at fil = 1.265(5) [30] to display the 
typical behavior already on quite small lattices. A sensitive measure of the strength of 
the transition is the order-order interface tension a„„ between the -I- and — phase, which 
turns out [29] to be (7„„ = 0.03443(47), 0.09785(60), and 0.16577(73) [31]. 

To evaluate the performance^ of iIk^ multicanoni("al nniltigrid algoritlim. w(^ r(x:or<led 
the time series for several ol)Sorval)los O and studietl their autocorrelation functions 
^U) — pU)/p(0)' where p{j) = {0,0,^j) — (O,)^. The exponential autocorrelation time 
^exp fQiiQ.5jrs from its asymptotic decay for large j, A{j) <x exp(— j/t"''''), and the integrated 
autocorrelation time r™' is obtained in the large k limit of T(fc) — 1/2 + Y^j=i A{j)- Here 
we shall concentrate only on the magnetization m which refiects most directly the (pseudo) 
dynamics of the transitions between the two ordered phases. In general we are interested 
in the variance of the (weakly biased) estimator O = Y,i'j!lw(mi)0,/ Y,^J"iw(mi) = 
ui,0,/Wi for a canonical expectation value (C')can- To facilitate a direct comparison with 
canonical simulations, we define for multicanonical simulations an (effective; autocorrela- 
tion time t"^ by the standard error formula for N„^ correlated measurements, 

= ^L^r^VN„,. (2) 

where a^^^ = (Of }caii — (Ci)fan is the variance of single measurements in the canonical 
distribution. The squared error can be estimated either directly by jack-knife block- 
ing procedures [32], or by applying standard error propagation to O = WiOi/Wi, which 
involves the (multicanonical) variances and covariances of WiOi and Wi, and the three 
associated integrated autocorrelation times r^.^ = Tq', t™^.^^ = t^q, and t™^.^ = fo-mO 
[28, 29]. 

Some results for /i^ = 1.30 from runs with at least 20 000 t"'^j Metropolis sw(M^ps or W- 
cycles (using piecewise constant interpolation without postsweeps) are collected iu Table 1. 
We see that r^^^ and r™' agree well with each other, indicating that the corresponding 
autocorrelation function A{j) can be approximated by a single exponential. While the 
agree with within error bars, the integrated autocorrelation times r^"^ are significantly 
smaller. This implies that for this observable^ A(j) is composexl of many diffen'ont modes. 
We further observe that the difference between r^'^^ and r^f can be quite appreciable, 
refiecting the varying probability distribution shapes with increasing L. For comparison 
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Figure 2: Effective autocorrelation times r^f for the model (1) in d = 2 with g — 0.25 as 
a function of lattice size L. 

with previous work we have also included the diffusion times r*'', which are defined as 
one quarter of the time it takes to travel from one (canonical) peak to the other and back. 

The log-log plot in Fig. 2 shows the behavior of r^f as a function of the lattice size 
L. Least-squares fits of to a pow(^r law. t''^ (x V" . yield for both update algorithms 
exponents of a 1.2,1.4, and 1.5 for — 1.30.1.35 and 1.40. i.e.. the multigrid update 
does not improve the asymptotic behavior. The off'octivo autoconolatiou times of tlie 
W-cycle, however, are about 20 times smaller than those of the Metropolis algorithm. 
If we finally take into account that a W-cycle requires more elementary operations than 
a Metropolis sweep, we obtain with our implementations on a CRAY Y-MP a real time 
improvement factor of about 10. 




Tabl(^ 1: Autocon'(^latioii tim(^s in units of sw(^(^ps r(^sp. (^'(^(^s for the Metropolis (M) and 
multigrid W-cycle (W) update in multicanonical simulations of the model (1) in d = 2 
with g = 0.25 and fj.^ = 1.30. 
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^exp 
^iiit 


212(12) 
204.4(4.0) 


11.30(32) 
10.88(12) 


668(23) 
690(11) 


37.2(2.0) 
34.69(76) 


3120(200) 
2984(63) 


148(11) 
150.0(4.0) 


746(62) 
758(37) 


_exp 
wm 
int 

wm 


209(12) 
171.1(3.4) 


11.34(33) 
9.82(11) 


655(31) 
509.8(8.9) 


36.9(2.0) 
27.58(59) 


2880(190) 
1840(40) 


146(13) 
96.6(2.4) 


600(120) 
374(23) 


' m 


322.7(6.1) 


18.51(20) 


1258(21) 


67.4(1.3) 


6050(120) 


321.9(7.6) 


1724(86) 


T-flip 

' m 


463.5(6.4) 


30.82(25) 


1759(24) 


91.7(1.3) 


7780(140) 


428.2(8.9) 


1922(85) 



6 



3 Refined Data Analysis 



3.1 Interface Tension 

A qiiandit.y of central importance for the kinetics of first-()r(l(^r plia,S(^ transitions is tli(^ 
interface tension a between the coexisting pliases [1. 2]. As discussed earlier, on finite 
periodic lattic(^s of size L' . this is r(^fi(M't(Ml by a doubl(^-p(^a,k stru(dur(^ of tli(! probability 
distribution for the energy or magnetization, with the minimum between the two peaks 
dominated by mixed phase configurations with two interfaces contributing an excess free 
energy of 2aL'^~^. This suggests [43] to extract the interface tension from the infinite 
volume limit of 

2'^*" = ;^M-Pi''x/^'i^]). (3) 

For accurate results the system lias to travel many Limes between the two peaks, which 
is a serious problem in canonical simulations of stron"' first-order transitions (large a). 
The multicanonical sampling, on the other hand, is just designed for this purpose since it 
gives the same relative errors for P^^], and P!^^l and thus optimizes the error on cr'^'. 

Using this so-called histogram method interface tensions have been estimated at the 
temperature driven first-order phase transitions of 2D g-state Potts models [9-12, 14, 51] 
and Nt — 2 SU(3) lattice gauge theory [10, 16], and at the field-driven transitions of the 
2D and 3D Ising [15] and the 2D (f>* model [29], For the 2D Ising model good agreement 
with the exact Onsager formula was obtained. For the 2D g-state Potts model with 
Hamiltonian [34] 

H^-jY.5,,,r, s. G{l,...,g}, (4) 

some uumerical results are compiled in Table 2 and compared with an analytical formula 
[35] (d(^riv(Ml iiftcr the first numerical results were already published), which follows by 
combining (a) tlie assumption of couiplete welting [36, 37] witli (b) duality arguments 
[38] and (c) an exact result for the correlation length id[lit) in the disordered phase at the 



Table 2: Comparison of analytical and numerical results for the order-disorder interface 
tension 2aod, in 2D g-state Potts models. 



Q 


4; 




2a„, (MC) 


7 


48.095907 


0.020792 


0.0241(10) .Janke (;/ al. [10] 
0.02348(38) Rummukainen [11] 
0.0228(24) Grossmann and Gupta [12] 


8 


23.878204 


0.041879 


0.045 Janke [51] 


10 


10.559519 


0.094701 


0.09781(75) Berg and Neuhaus [9] 
0.10 Janke [51] 
0.0950(5) Billoire et al. [14] 


15 


4.180954 


0.239179 


0.203(9) Gupta [42] 


20 


2.695502 


0.370988 


0.3714(13) Billoire et al. [14] 
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transition point fit ~ ln(l -|- ^) [39] (see also [40, 41]), 



2a^ = o„„ = lli, = -f:\n\\±^\, (5) 

where w„ = [v^cosh ((/) + i)7rV2ti)]"^ with v = In (i [y/^-b 2 + 2j). More 

precisely, it was shown [3G] that 2acd < f^oo for «// q > 5. The opposite ine<iuality could 
only be d(Tiv(Ml for q > q^ (with 4 < qg < OC' b(4ng a sut(ici(nitly larg(^ constant), but by 
basic thermodynamic arguments it is commonly believed that it actually also holds for 
all q>5. 

Overall the numerical and analytical values are in good agreement, but noteworthy is 
the systematic trend of the numerical data to overestimate the analytical values, which 
are actually exact upper bounds. 

3.2 Correlation Length 

To test the formula (5) for the correlation length (d(pt) directly, we performcMl further 
Monte Carlo simulations [44] of the Potts model (4) and measured the ky = projection 
g of the correlation function 

G{tJ) = {S.,., - 1/g), (6) 

at l3t in the disordered phase. Here we simply choose the lattice sizes large enough in order 
to suppress transitions to the ordered phase. Estimates of autocorrelation times showed 
[44] that in this situation the optimal update procedure consists of many single-cluster [45] 
update steps combined with one multiple-cluster [46] update for efficient measurements 
using the improved estimator 

G{i,j)^^—^{G{^,j)), (7) 

where &{i,j) = l, Hi and j belong to the same cluster, and Q = otherwise. 

The numerical data for q = 10 al /?f are shown in Fig. 3. From the curvature of g{x) 
for small distances it is clear that fits with the simplest Ansatz for periodic boundary 
conditions. <j{:r) — acosh((:r — L/2)/(d), will only work for very large x. We therefore 
tried to use the more general Ansatz ^(a;) = a cosh((2; — L/2)/^d) + 6cosh(c(2; — L/2)^d), 
with four parameters a, b, c, and ^d- The solid line in Fig. 3 shows a fit where ^d is held 
fixed at its theoretical value and only the three remaining parameters are optimized. The 
good quality of the fit over a wide x range shows that our data are compatible with 
the tli(X)r(4.ical (^xp(M'lations. To get a more quantitative measure we next performed 
four-parameter tits witli as a free parameter. Fits over the same x range yield an 
estimate of (,;(/it) = 9.2(8), which is about 15% smaller tliau tlie exact value but within 
the statistical (^rrors still consistent. By restricting the fit int(^rval to larg(^r :r values, we 
observed a tendency to higher values but with the drawback of increased error bars. Since 
our data for q = lb and 20 exhibit the same qualitative behavior, we believe that higher 
excitations cannot be neglected at the distances studied so far (x^ax = L/2 Ri 7^d)- To 
cope with this problem we are currently performing further siniTdations on 2L x L lattices, 
which should allow to study the correlations over even larger distances. 
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Figure 3: Semi-lo"' plot of the correlati(m function fj(x) vs distance a; for g = 10. The 
solid line is a three-parameter fit with held hxed at its theoretical value. 

3.3 Accurate Determination of the Transition Point and La- 
tent Heat 

The traditional way to locate first-order transition points is based on the scaling boliav- 
ior of the speciflc-heat maxima C^i^^{V) and Binder-parameter minima Bmiii(^'')- wL(T(^ 
CiVJj) ^ fPV({e^)-{e)^) and B(V,l3) = I- {e*) . In the idealized infinite volume 

limit. tli(V(^ quantities develop singularities at the transition point, In a finite volume 
V = L"^ the singularities are smoothed out and, depending on its strength, the transition 
is signalized by more or less pronounced finite peaks of the specific heat or dips of the 
Binder parameter near the infinite volume transition point. For an illustration see Fig. 4 
where data for the two-dimensional 8-state Potts model are shown. If the volume is cubic 
or n(^arly cubic tli(^ location of tli(^ (^xlnnna is lypically sliift(Hl by a,n amount (){V~^) 
with respect to the actual infinite volume transitiou point and one may try to estimate 
fit = l/Tt from the finite volume results by extrapolations in 1/V [47]. 

In addition to random statistical errors the data are in general also systematically 
shifted by exponential corrections which are difficult to take into account theoretically. 
This makes the extrapolation of finite volume data not always reliable. It is therefore 
desirable to find definitions of a finite volume transition point which involve no power-law 
corrections at all. 

The starting point for such dotiuitious is the observation that on lattices with periodic 
boundary conditions the partition function of a, model describing the co(^xist(nice of one 
disordered and q ordered phases can be written for large enough q as [48, 49] 

Z,.AV, /3) = f E e-»-('^>A (l + O (Fe-^/^")) . (8) 

Vm=0 / 

Here Lq < °o is a constant, L is the linear length of the lattice, and f„i{p) is the meta- 
stable free energy density of the phase m. It can be defined in such a way that it is equal 
to the idealized infinite volume free energy density /(/^) if m is stable and strictly larger 
than if m is unstable [48, 49]. This implies (see also Refs.[50, 51]) that in the infinite 
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volume limit the parameter 

Ar(y,/3) = Zp„(y,/3)e«<«^ (9) 

is equal to the number of stable phases at the inverse temperature /}, i.e., N{l3) = 

limy f^N{V,f)) = q in the ordered phase, N{f)) = q + 1 a.t the transition point /3t, 

and N{l3) = 1 in the disordered phase. A natural definition of a finite volume transition 
point /3t(V) is thus the point where N{V,j3) is maximal. Due to the bound (8) (and 
similar bounds for derivatives [48, 49]) this definition leads to only exponentially small 
shifts of lit(y) with respect to the infinite volume transition point /?,. 

The free energy /(/j) in (9) is only defined in the tliormodynamic limit and hence 
not acc(vsibl(^ to iium(n'ical simulations. It is tli(n'(4br(^ ii(M'(^ssai'y to (4imiiiat(^ this U^xm 
by, e.g., forming a suitable ratio. Instead of (9) one may look for the maximum of the 
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number-of-phases parameter 



(10) 

where a = V2/V1 [50, 51]. By inserting (8) it is straightforward to verify that with 
increasing temperature N {Vi,V2, l3) smoothly interpolates between the values q, q + 1, and 
1. Actual simulation results for a = V2/V1 1.6 are shown in Fig. 4. The locations of 
the maxima define the desired, only exponentially shifted finite-volume transition points 
/3((T-'i-l''2) which will ho doiioted l3v/v- By differentiating InA^ with respect to /J one 
readily S(^(^s that <l(M,(n'miuiug l3v/v amounts to solving aE{Vi, f^v/v) = E(V2,l3viv) or 
e{Vi,l3v/v) = <'-{V2, l3v/v), i-e., to locating the crossing point of the internal energies per 
site, e = E/V, of the two lattices of different size. 

The numerical determination of f)v/v requires simulations on two different lattices. In 
Ref.[50] we have therefore proposed another definition of a finite- volume transition point 
which requires data from one lattice only. Its definition is based on the fact that at the 
infinite-volume transition point all free energies are equal, so that eq. (8) implies 

w„=J2 e-*^"'*'^ = qe-^'f'^f'^^ = qw,, (11) 

m— 1 

apart from exponentially small corKM'tions. Here the free energy density of the "zeroth", 
disordered phase is douotod by /j, and it!„^j are the associated statistical weights of the 
coexisting phases. A natural definition of a finite-volume transition point is thus the 
point where the ratio of the total weight of the q ordered phases to the weight of the dis- 
ordered phase approaches q. More precisely we introduce the ratio-of-weights parameter 

R(V,f)) = WJW,= J2 PvAE)l E PvAE), (12) 

where PvA^) are the (double-peaked) energy histograms, and determine fiw by solving 

R(V,M = q- (13) 

The parameter Ea iu (12) is defined by rowoiglitiug [52] the energy distribution to the 
temperature where the two p(^aks of Pvji(E) liav(^ (Hjiial h(in'lit and then taking Eq as 
the energy at the minimum between the two peaks. Other dotinitions of Eq would be 
reasonable as w(41. a,s for example tli(^ int(n'ual (^k^'^'v a,t th(^ t(^nperatiir<! where? the specific 
heat is maximal. Since it is expected that the relative height of the minimum between 
the two peaks decreases like exp(— 2(7L'*^^) as L — > 00, all these definitions do in fact 
only differ by exponentially small corrections and it is a matter of practical convenience 
to choose Eq. In Fig. 4 the logarithm of the ratio-of-weights parameter R(V,li) is plotted 
as a fun("tioii of t(inperature. 

In (13) wo liavo asRiimed that the iiunihor of ordered phases, </, is known by general 
arn'iiiiKnits. If this is not the cas(\ oik^ may iis(; the crossing points jiw/w SH.tisfyiug 
R {yi,Pwlw) = R (T^2,/^w/w) as estimates for pQ. The value of R at the crossing point 
then gives the ratio of the number of coexisting ordered and disordered phases. This, 
however, requires again the simulation on two lattices of different size. 



N(VuV2,l3)-- 



Z,AV2ji) 
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(a) 



(b) 



0.000 0.001 0.002 0.003 0.000 0.001 0.002 0.003 0.004 0.005 

1/V 1/V 

Figure 5: Scaling behavior of the finite-volume transition points dotinod by the Binder- 
parameter minimum (o), the specific-heat maximum (•). tli(^ maximum of the number- 
of-phases parameter (A), and the ratio-of-woiglu. parameter (□). Tlio solid straight, linos 
are the exactly known IjV corrections corresponding to o and •, and the dashed, almost 
interpolating curves show exponential fits (including the IjV corrections) to these data. 
Note that the (1/y)^ corrections are almost invisible on this scale and in any case point 
in the "wrong" upward direction. The long dashed horizontal lines indicate the exact 
infinite-volume transition points. 



The convergence properties of the four finite-volume transition points are compared 
in Fig. 5 for a very weak (q = 5) and a rather strong (q = 8) first-order transition. Notice 
the surprisingly fast convergence of l3w{y) even for g = 5. 

The ratio-of-w(nghts method leads naturally to a finit(^-volume definition of the latent 
heat [51] which also should have only exponentially small corrections with respect to the 
infinite volume limit. Since 



ln(»„/™,) = -/3y (/„-/,), (14) 
the slopes of R(V,f)) in Fig. 4 at the crossing point may be used to define 

Ae{V) = e,(V)- e,{V) = ln{WJW,)/V = HWJW,)/V. (15) 

The resulting estimates Ae{V) are plotted in Fig. 6 and compared with the traditional 
definition based on the peak locations of PvA^) [53]. For strong first-order transitions 
((? = 8 and 10) the asymptotic limit is indeed reached much faster with the new definition. 
For a very weak transition (g = 5), on the other hand, both methods yield comparable 
estimates which are still far away from the limiting value. 
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Figure 6: The finite-volume lat(^it L(^at Ac of tL(^ 2D q-Mr,iic Pol.ts mo(l(4 vs liiK^ar lal.tic(^ 
size L. The open symbols show the t.radit.ioual estimates from the peak locations of 
Pv,i}(E), and the filled symbols follow from tlio slopes of the ratio-of- weight parameter. 
The dashed horizontal lines show the exactly known infinite- volume limits [54], 



4 Concluding Remarks 

TL(^ mat(n'ial r(^vi(nv(Ml in these lecture woUv is. a(lmitt(Mlly. a sonKnvhat biased selection 
of recent developments in this field. Different improved Monte Carlo update schemes 
are discussed, e.g., in Refs.[55, 56]. Also, refined data analyses have been developed 
along many different lines. For a recent review focusing on the analysis of probability 
distributions and various cumulants, e.g., see Ref.[57]. 

The combination of multicanonical reweighting with multigrid update schemes leads 
to a significant reduction of autocorrelation times by a roughly constant factor. It still 
remains a challenge to develop Monte Carlo algorithms which exhibit no slowing down at 
a first-order phase transition. 
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